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ABSTRACT 

Results of a simulation of synchrotron-self Compton (SSC) emission from a 
rotation-powered pulsar are presented. The radiating particles are assumed to 
be both accelerated primary electrons and a spectrum of electron-positron pairs 
produced in cascades near the polar cap. They follow trajectories in a slot gap 
using 3D force-free magnetic held geometry, gaining pitch angles through reso¬ 
nant cyclotron absorption of radio photons, radiating and scattering synchrotron 
emission at high altitudes out to and beyond the light cylinder. Full angular 
dependence of the synchrotron photon density is simulated in the scattering and 
all processes are treated in the inertial observer frame. Spectra for the Crab and 
Vela pulsars as well as two energetic millisecond pulsars, B1821-24 and B1937-I-21 
are simulated using this model. The simulation of the Crab pulsar radiation can 
reproduce both the hux level and the shape of the observed optical to hard X-ray 
emission assuming a pair multiplicity of M_|_ = 3 x 10^, as well as the very-high- 
energy emission above 50 GeV detected by MAGIC and VERITAS, with both 
the synchrotron and SSC components reflecting the shape of the pair spectrum. 
Simulations of Vela, B1821—24 and B1937-I-21, for M+ up to 10^, do not produce 
pair SSC emission that is detectable by current telescopes, indicating that only 
Crab-like pulsars produce signihcant SSC components. The pair synchrotron 
emission matches the observed X-ray spectrum of the millisecond pulsars and 
the predicted peak of this emission at 1 - 10 MeV would be detectable with 
planned Compton telescopes. 


1. Introduction 

The recent detection of pulsed emission from the Crab and Vela pulsars at energies 
above 50 GeV by ground-based air-Cherenkov telescopes marks a turning point both in 
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technological achievement of the detectors and in pulsar acceleration and radiation modeling. 
The MAGIC telescope hrst detected pulsed emission from the Crab above 25 GeV (Alin et 
ah 2008), followed by pulsed detections above 100 GeV by VERITAS (Alin et ah 2011) and 
up to 400 GeV with MAGIC (Aleksic et ah 2012). The averaged and phase-resolved spectra 
above 25 GeV are consistent with a smooth broken power-law extension of the spectrum 
measured by the Fermi Large Area Telescope (LAT) above 100 MeV, thus ruling out an 
exponentially cutoff power-law spectrum predicted by curvature radiation (CR) models for 
GeV emission (e.g. Romani 1996, Harding et ah 2008 [HSDF08]). Very recently, pulsed 
emission was detected from the Vela pulsar above 30 GeV with the HESS (Stegmann et ah 
2014) telescope and above 50 GeV with the Fermi-LAT (Leung et ah 2014). In this case 
it is not clear whether another emission component in addition to CR is required. Most 
recently, VERITAS has reported a non-detection of the Geminga pulsar above 100 GeV 
(Alin et ah 2015). McCann (2015) performed a stacking analysis of 115 Fermi pulsars from 
the Fermz-LAT 2nd Pulsar Catalog (Abdo et al. 2013), both young and MSPs, and found 
no signihcant emission at more than 7% Crab level above 50 GeV. 

For many years, CR was the favored mechanism for the high energy (> 100 MeV) pulsed 
emission in a number of different models for pulsar emission, and VHE emission similar to 
what has been detected from the Crab was never predicted. In models where acceleration 
and emission occur above the pulsar polar caps (PC models, Daugherty & Harding 1996), 
the predicted spectrum of CR from primary electrons has a very sharp “super-exponential” 
cutoff due to attenuation by the magnetic pair production process that occurs in the very 
strong magnetic fields near the neutron star surface. Measurement of the spectral cutoff of 
the Vela pulsar by Fermi (Abdo et al. 2009) has shown that a super-exponential cutoff can 
be eliminated with high significance, ruling out emission from near the neutron star PCs. 
Models for emission originating in the outer magnetosphere are therefore favored, since many 
of these predict curvature radiation spectra from primary particles with exponential spectral 
cutoffs at energies of a few GeV that agree well with those measured for most pulsars by 
Fermi. Such models include the outer gap (OG), where acceleration and emission occur 
in vacuum gaps that form above the null charge surface (e.g. Romani 1966, Cheng et al. 
2000, Hirotani 2006) and the slot gap (SG), where acceleration and emission occur in a 
narrow region along the last open field lines from the neutron star surface to near the light 
cylinder (Muslimov & Harding 2004 [MH04], HSDF08). More recently, models postulating 
high-energy emission originating outside the light cylinder, near or in the current sheet that 
forms in the spin equatorial plane have been developed. Several of these models assume 
that synchrotron radiation (SR) from particles accelerated by magnetic reconnection in the 
current sheet produces the GeV emission (e.g. Petri 2012, Uzdensky & Spitkovsky 2014), 
while others assume that CR by particles accelerated by induced electric helds in dissipative 
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magnetospheres produces the GeV emission (Kalapotharakos et ah 2014, Brambilla et ah 
2015). 

The Fermi LAT has to date detected more than 160 7 -ray pulsars and all of these with 
statistics adequate for spectral htting of a power law with exponential cutoff 

dN{E)/dE = E-^ exp[-(E/Ee)^], (1) 

show cutoffs, Ec, in the range 1-7 GeV (Abdo et ah 2013), assuming 6 = 1 . For a subset of 
pulsars having higher counts, hts to the above expression with b free can be performed and 
a number of these are best £t with 6 < 1 , including the Grab and Vela pulsars, indicating 
a more gradual cutoff than pure exponential. This finding raises the question of whether 
the more gradual high-energy cutoffs are due to a blending of GR spectra from a number of 
different rotation phases (Abdo et al. 2010) and/or locations in the magnetosphere (Leung et 
al 2014), or additional emission components. It has been suggested that the very-high-energy 
(VHE) emission from the Grab is inverse-Gompton emission (Lyutikov et al. 2012), since it 
is difficult to produce photons above 100 GeV with GR. Lyutikov (2013) modeled the Grab 
X-ray and 7 -ray spectra as cyclotron-self Gompton emission from electron-positron pairs in 
the outer magnetosphere, fitting the data to determine the spectrum and multiplicity of 
the pair spectrum as well as the location of emission. An synchrotron-self Gompton (SSG) 
model for the Grab VHE emission from the outer gap (OG) was suggested by Aleksic et 
al. (2011) and Du et al. ( 2012 ) modeled the observed Grab emission from the annular gap. 
Alternatively, it was proposed that particles accelerated by reconnection in the current sheet 
reach temperatures of 10 GeV, if equipartition is assumed, and their SR can extend to 100 
GeV in the observer frame after Doppler boosting (Uzdensky & Spitkosky 2014, Iwona & 
Petri 2015). 

In this paper, we apply a pair cascade simulation coupled to an outer magnetosphere 
radiation model to predict the broadband X-ray to VHE 7 -ray spectra of several classes of 
pulsar. In this model, the pairs are produced in cascades above the PGs initiated by primary 
accelerated electrons and sustained by magnetic pair creation (e.g. Daugherty & Harding 
1982). The pairs lose their momentum perpendicular to the magnetic held to SR near the 
PGs and then stream into the outer magnetosphere, where they absorb radio photons in the 
cyclotron resonance to regain pitch angles and emit SR (HSDF08). The SR provides soft 
photons for inverse-Gompton scattering by these same particles to boost their energies to the 
7 -ray range. We also include GR radiation of primaries accelerated by electric fields in a SG 
geometry, as well as their SR components due to cyclotron resonant absorption. The radio 
emission is simulated using an empirical cone and core beam model. Since the magnetic held 
structure in the outer magnetosphere becomes important for these models, we use a force-free 
magnetosphere configuration to define the magnetic held, a significant improvement over the 
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retarded vacuum dipole field that we had used previously (HSDF08). Use of the force-free 
magnetic field structure allows us to follow particles and radiation arbitrarily close to and 
even beyond the light cylinder, which should not be a physical limitation of emission models. 
Although force-free magnetospheres assume that the electric field parallel to the magnetic 
field, U||, vanishes, they well approximate the field structure of energetic pulsars where 
departures from force-free conditions exist only in narrow accelerator gaps. We therefore 
use the force-free field structure but assume acceleration of primary particles in a narrow 
SG. Using these simulations, we can address the question of how many pulsars should have 
detectable SSC emission and how high in energy this emission extends. 


2. 3D Magnetosphere Geometry 

Since radiation of relativistic particles is beamed along their direction of momentum, and 
particle trajectories in pulsar magnetospheres follow closely (but not exactly) the magnetic 
field lines, the structure of the magnetic field is a critical component to emission models. A 
static dipole is not a good approximation to a pulsar magnetic field, even near the neutron 
star, since the rotation causes a sweepback of the field lines near and beyond the light 
cylinder. The vacuum retarded dipole solution (Deutsch 1955), adopted in many radiation 
models (e.g. Romani & Yadigaroglu 1995, Cheng et ah 2000, HSDF08), incorporates the 
sweepback due to retardation of the field and produces a distortion and shift of the open 
field volume (including the PC) toward its trailing edge (Dyks & Harding 2004). Numerical 
solutions of force-free magnetospheres (Spitkovsky 2006, Timokhin 2006), where the plasma 
density is assumed to be high enough to screen E\\ so that the ideal MHD condition E ■ B = 0 
is enforced, show a larger degree of magnetic field sweepback, as well as straightening of the 
poloidal field lines due to currents. Recent simulations of dissipative pulsar magnetospheres 
(Kalapotharakos et ah 2012, Li et al. 2012), that drop the ideal MHD condition and allow 
finite E||, show that the structure of the magnetic field stays close to force-free for the high 
plasma conductivity needed to match the y-ray light curves (Kalapotharakos et al. 2014, 
hereafter KHK14) and spectra (Brambilla et al. 2015) of young pulsars. 

We will therefore adopt the force-free magnetosphere structure to model the particle tra¬ 
jectories and radiation, and all calculations are done in the non-rotating, inertial observer’s 
frame (lOF). The pairs are assumed to be produced at the PCs and to screen the E\\ at 
higher altitudes, except for a narrow gap along the last open field lines (the SG). They will 
therefore not experience any acceleration in the calculation of their radiation in the global 
magnetosphere. The primary electrons are assumed to be accelerated only in the SG, with a 
constant E\\ that we take as a free parameter. Both pairs and primary electrons are injected 



5 


at the neutron star surface at footprints of open field lines that are determined by “open vol¬ 
ume coordinates” (Dyks et ah 2004 [DHR04]). These coordinates map the open held region 
bounded by held lines that close within the light cylinder, i?LC = c/hl, where hi is the pulsar 
angular rotation frequency. The magnetic held vector is determined at each point by 3D 
interpolation of a table of values read-in from a numerical force-free magnetosphere solution 
(Kalapotharakos et ah 2012). This solution has a resolution of 0.02i?Lc and extends to a 
minimum radius of 0.2i?LC) which is the surface of a pulsar with period 0.8 ms. To extend 
this solution to the surface of pulsars with longer periods, we join the numerical solution to 
a retarded vacuum dipole solution over the range 0.2 — 0.4i?LC using a ramp function. 

The code hrst computes the PC rim at the neutron star surface by performing 4th order 
Runge-Kutta integrations along the held lines of the numerical solution to determine the open 
held footpoints. This is done by iteration, choosing an initial value of magnetic polar angle 
di at a number of diherent values of magnetic azimuth 0. If the held line with that footpoint 
does/does not close within Rlc, a smaller/larger value of 6 is chosen until the closure point is 
at Rlc within a given tolerance. We dehne “open volume coordinates” (rovc^ovc) to identify 
footpoints at the neutron star surface of held lines along which the particles are traced. The 
“radial” coordinate rove, equivalent to magnetic polar angle, is equal to 0 at the magnetic 
pole and 1 at the PC rim. Lines of constant rove dehne a set of concentric deformed rings 
that hll the PC surface (see Figure 2 of DHR04). The “azimuthal” coordinate /qvc measures 
the arclength along each distorted ring with hxed rove- The /qvc increase in the direction of 
magnetic azimuth, c/pc, iu a counterclockwise direction around the PC, starting from /qvc = 0 
at (()pc = 0, dehned to be at the magnetic meridian (i.e. the line between the magnetic and 
spin axes). With this dehnition, (ppc = 37r/2 is on the leading side of the PC and (ppc = n/2 
is on the trailing side (Note that (ppc = 0 -f tt). 

To model an injection of particles that is uniform over the PC, we use a set of rings 
between r™™ and with equal spacing dove- These footpoints are then spaced uniformly in 
each of the rings with Ni = Aazim360 equal divisions. Since each ring has the same number of 
footpoints, and the rings have varying circumferences, the contribution from different rings 
must be weighted by /ring/4im,where /nng is the length of a particular ring and l^ra is the 
total PC rim length. To determine the trajectories of particles in the lOF, we require the 
total particle velocity to be the sum of a drift component and a component parallel to the 
magnetic held (KHK14), 

f ExB , 


V = 


( 2 ) 
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where (Gruzinov 2008, Li et al. 2012), 


El = 

( 3 ) 

Bl = (B^ - + v/(B2-E2)2 + 4(B-E)2)/2 

( 4 ) 

Bo = sign(B ■ E)JbI. 

( 5 ) 


By requiring that v ~ c and that the motion is outward, we can solve for the scalar quantity 
/ at each point along the trajectory. In the special case of a force-free magnetosphere, v 
stays parallel to B in the corotating frame. In the lOF, the step size along the particle 
trajectory is M = ds /where ds is the step size along the rotating magnetic field line (see 
Appendix C of Bai & Spitkovsky 2010). In treating the particle trajectories in the lOF with 
force-free held geometry, it is possible to extend the trajectories and radiation outside the 
light cylinder. Although the particles are following the held, they slide forward along the 
held lines as the held lines sweep back, following a nearly radial path in the lOF at r > i?LC- 


3. Pair Cascade Spectra 

The process of electron-positron pair production in pulsar magnetospheres is thought 
to be critical for generating charges for the magnetosphere and for the pulsar wind nebula, 
as well as plasma necessary for coherent radio emission. The pairs are produced in elec¬ 
tromagnetic cascades above the PCs (Daugherty & Harding 1982) or in OGs (Cheng et al. 
1986). Here, we assume that the pairs that radiate SR and SSC in the outer magnetosphere 
originate from PC cascades, that are initiated in the strong electric helds near the neutron 
star surface by acceleration primary electrons. We have calculated the spectra of pairs using 
a code that simulates a steady (non-time-dependent) electromagnetic cascade above the pul¬ 
sar PC (Harding & Muslimov 2011 [HMll]). Primary particles are accelerated by an electric 
held induced by rotation of the magnetic held, derived assuming space-charge-limited how 
(i.e., free emission of particles from the neutron star surface) (Arons & Scharlemann 1979), 
and emit curvature radiation. The highest energy curvature photons are absorbed by mag¬ 
netic pair attenuation (Erber 1966, Daugherty & Harding 1983), producing a spectrum of 
hrst-generation electron-positron pairs. The pairs are born in excited Landau states (with 
non-zero pitch angles) and radiate synchrotron photons that spawn further generations of 
pairs. Since the total cascade multiplicity M+ (average number of pairs produced by each 
primary particle) depends on the pulsar period P and surface magnetic held strength Rg, 
younger pulsars produce high M+ cascades but older pulsars, with low magnetic helds and 
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long periods, have much lower M_|_. This results in a pair death line in P-P which, in the 
case of a dipole held, predicts that the older half of the pulsar population cannot produce 
signihcant pair multiplicity. 


However, the sweepback of magnetic held lines near the light cylinder (due to real and 
displacement currents), as well as possible held distortions near the neutron star will cause 
the magnetic PCs to be ohset from the dipole axis (e.g. Dyks Sz Harding 2004). HMll 
introduced a perturbed dipole held structure with non-axisymmetry and showed that the 
resulting ohset PC enhances pair multiplicity by increasing the accelerating electric held. 
For the pair cascade simulations, we will adopt this magnetic held structure which has two 
conhgurations for the dipole ohset in which the magnetic held is symmetric or asymmetric 
with respect to the dipole axis. Modeling of the thermal X-ray light curves of several MSPs 
requires an ohset of the PC hotspot that is symmetric (Bogdanov et al 2013). We adopt the 
symmetric distorted held in simulating the pair spectra of the MSPs in this paper. In this 
case, the magnetic held in spherical polar coordinates {rj, 6, 0 ) is 


B — 


f cos^ + - 6 {1 + a) sin^ — 0 e sin 6 ' cos 6 ' sin((;/) — ^o) 


( 6 ) 


where r] = r/R is & dimensionless radial coordinate in units of the stellar radius R, Bg is the 
surface held strength at the magnetic pole, a = e cos (0 — 0 o) is a parameter dehning the 
azimuthal distortion of the polar held lines from a dipole, and 0 o is the magnetic azimuthal 
angle that dehnes the plane of the PC ohset. HMll derive the parallel component of the 
electric held, E\\, using this held structure and we have used the of Eqn (11) of HMll 
that corresponds to a symmetric ohset. We use this E\\ to accelerate the electrons from 
diherent starting points over the PC surface and simulate the pair cascades over the whole 
PCfor a range of P, P (or equivalently, Bg), and ohset parameter e. We assume a non-zero 
ohset only for the MSPs. Pair spectra for the parameters of the Crab (e = 0) , Vela (e = 0) 
and MSPs PSR B1821-24 (e = 0.6, (/)o = 37 r/ 2 ) PSR B1937+21 (e = 0.6, 0o = 37 r/ 2 ) are 
shown in Figure The Crab pair spectrum extends over hve decades of pair energy, from 
7 _i_ = 20 to 7 + = 10®, while the Vela pair spectrum has a smaller range, cutting oh around 
7 _i_ = 4 X 10®. The diherence is due to the shorter period and larger PC size of the Crab, 
which gives smaller radii of curvature of the held lines, so that CR photons have higher 
energy and gain large angles to the held, and the higher magnetic held strength. The pair 
spectra of the MSPs start at energies about two decade higher, 7 + ~ 2 x 10® and extend to 
higher energy, 7 _|_ ~ 10^. This striking diherence in pair spectra of young and MSPs comes 
from their large diherence (four orders of magnitude) in surface held strengths. As described 
by HMll, the photons must have higher energies to produce pairs in lower magnetic helds, 
which raises the minimum pair energy, but are also able to escape at higher energies, raising 
the high energy cutoh in the pair spectrum. The non-zero PC ohset decreases the lower 




cutoff of the pair spectrum relative to the pair spectrum with zero offset by a fraction of a 
decade. 

Simulations of PC cascades including time dependence of the electromagnetic fields 
and plasma production (Timokhin 2010, Timokhin & Arons 2013) hnd that to satisfy the 
charge and current density of the global magnetosphere models the pair cascades are non¬ 
steady. The accelerating electric held is time-dependent, in the general case where the 
magnetospheric current density J is not very close to the Goldreich-Julian value, Jgj = Pgjc, 
with cycles of pair creation followed by complete screening of the electric held. The timescale 
of these pair cascade cycles is very short (of order the light crossing time across the gap). 
Snch non-steady pair cascades can produce pair multiplicities that are much higher than 
the steady cascades. For young pulsars such as the Crab and Vela, the time-averaged pair 
multiplicities are about several times 10^ (Timokhin & Harding 2015), compared to ~ 2 x 10^ 
for the Crab and ~ 6 x 10^ for Vela steady pair cascades. Time-dependent pair cascade results 
are not yet available for MSPs, where 2D simulations are reqnired. Althongh it is beyond 
the scope of this paper to simulate time-dependent pair spectra (especially for millisecond 
pulsars [MSPs]), we will allow for pair multiplicities higher than that of steady cascades in 
our SSC simulations. 


4. Simulation of Emission 

We model the high-energy radiation over the entire spectrum from optical to VHE 7 - 
ray wavelengths. Emission from both primary particles and pairs are simulated as they 
move along their trajectories with step size A£, starting at the neutron star snrface to a 
maximum radius rmax- Primary electrons are assumed to undergo continuous acceleration 
with a constant electric held, d'y/di = R^cc in a narrow gap along the boundary of the open 
held region, between Vovc = 0.95 — 0.99. The electric held in the high-altitnde gap is not the 
same held from HMll that is used for the pair cascades, since the HMll held is only valid for 
low altitudes near the PC. The high-altitude gap forms at the outer edge of the PC, where 
the interior Ey decreases to valnes too low for pair cascades to develop. Although analytic 
solntions exist for both the low altitude SG (Muslimov & Harding 2003) and its extension to 
high altitudes (MH04), we have chosen not the use the solution of MH04 for several reasons. 
First, the gamma-ray luminosity modeled using the SG solution of MH04 with widths 
narrow enough to produce the observed light curves falls short of the observed luminosities 
by abont a factor of 10 for young pnlsars. As noted by Pierbattista et al. ( 2012 ), the Ey 
solntion for the original SG and also the OG do not provide enough gamma-ray luminosity 
to match the observed Fermi pulsars. The Eh of HM2011 on the offset side of the PC is 
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larger and may produce a large enough at high altitude to produce luminosities matching 
the data. However, it is very difficult to derive the E\\ extension for the offset PCs and 
this has not yet been done. Also, with the advent of dissipative magnetosphere models 
we will soon be able to have more self-consistent E\\ throughout the magnetosphere, albeit 
numerical. Until these are available, we decided it is better to treat E\\ at high altitudes as a 
free parameter. We have therefore simply assumed a constant E\\ for the high-altitude gap, 
similar to the behavior of the high-altitude symmetric SG MH04. 


The electron-positron pairs are assumed to experience no acceleration as they flow along 
their trajectories in the screened region inside the gap, between Tovc = 0.91 — 0.95. CR, SR 
and SSC radiation are simulated for all particles in the same way, as described below. The 
primary electrons are injected at the neutron star surface with very low Lorentz factors, 
7 = 2 , while the pairs are injected with the pair cascade spectra shown in Figure [Tj All 
particles are assumed to initially have momentum only parallel to the magnetic held, since 
the pairs lose all their perpendicular momentum to SR very near the neutron star surface in 
the pair cascade. Single primary electrons are injected at the start of each primary trajectory 
at the neutron star surface and their radiation is normalized to the primary hux 


h, = ncj cttR^ elc [{rZ:f - (7) 

where rioj = i?o^/2vrec is the Goldreich-Julian density and 6pc = is the PG 

half-angle. The primary particles from the interior of the PG are on held lines that do not 
thread the high-altitude accelerator gap, but are on held lines where the accelerating held is 
screened above a pair formation front by pair cascades at low altitude. Since these primaries 
are not accelerated above the pair front, they lose most of their energy quickly and do not 
contribute much to emission at high altitudes. The emission from the PG pair cascades is 
all emitted at low altitudes, below 6-7 stellar radii, and will only be visible when an observer 
viewing angle is very close to the magnetic axis. Since most observer angles will not see the 
PG emission, we therefore can neglect their contribution in the present study. 


To simulate radiation from a spectrum of pair energies, shown in Figure [T| the pair 
spectrum is divided into logarithmic energy intervals, A 7 +, between and 7 ™“* and a 
single pair member (electron or positron) at each energy, 7+, is injected at the base of each 
pair trajectory. Their radiation is normalized to the hux of pairs in each energy interval. 


’ 7 +( 7 +) 


2M+n+(7+)A7+hp 


^max 

n+(7+)d7+ 


( 8 ) 


and M+ is the pair multiplicity. We have kept M+ as an independent quantity rather than 
derive it from the integral of the pair spectrum so that it can be treated as a free parameter 
in the simulation. 
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The direction of photon emission, rjem, is along (3 = \/c, the direction of particle motion 
in the lOF, with photon emission angles 


f^em 


13 z, 


(pem — cLtcLll 



( 9 ) 


As the particle moves a distance di, the change in phase is A0rot = ^di/c = Qdt. The 
emitted photons are accumulated in sky maps of observer angle Cobs = acos(/iem) versus 
phase 0obs; both with respect to the pulsar spin axis (the z axis on the lOF Cartesian grid). 
To determine the observed phase, we apply the rotation and add the time delay correction, 

0obs = 0em - A^rot + (10) 

-Klc 

where rgm is the radius of photon emission. The quantities in the sky maps are the emitted 
photon fluxes, normalized using Eqn ([^ for primaries and Eqn ([^ for pairs, divided by the 
solid angle of that sky map bin, Af2 = sin Cobs ACobsA0obs- The phase-averaged observed flux 
at a viewing angle Cobs is then obtained by summing the fluxes in 0obs and dividing by 2TTd‘^, 
where d is the distance to the source. 


4.1. Curvature Radiation 


The energy spectrum of curvature radiation from a single electron with Lorentz factor 


7 IS 


Ncr{£) = \/3 — 7 k f^ 

C \^cr J 


where e is the emitted photon energy in units of mc^ and 

3 c o 

Scr = -—'y , 

^ pc 

and the function hi{x) is dehned as 

poo 

k{x) = X K^i2,{x')dx . 


( 11 ) 


( 12 ) 


(13) 


The radius of curvature pc in the force-free magnetosphere is not that for a pure dipole held, 
but is the radius of curvature of the particle trajectory determined in the lOF by computing 
the inverse of the trajectory curvature using the second derivative at the particle position: 


d& 


P. 


-1 


(14) 
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The approximate form of the photon spectrum is a power law with an exponential cutoff 


at e 


crt 


Ncr{£) 

de 


a 

(Amc)h3 


/ \ 


(15) 


4.2. Synchrotron Radiation 


Since all particles start their trajectories with momentum parallel to the magnetic held, 
they will not radiate any SR until they acquire hnite pitch angles. We assume that particles 
gain pitch angles by resonant absorption of radio photons (Shklovsky 1970, Lyubarski & 
Petrova 1998, LP98). In this process, relativistic particles absorb photons that have energies 
at the cyclotron resonant frequency in their rest frame, which results in an excitation of the 
particle to a higher Landau state and an increase in pitch angle. The particle will then emit 
spontaneous cyclotron radiation if its momentum perpendicular to the magnetic held is non- 
relativistic in the frame where the parallel momentum vanishes, or synchrotron radiation if 
the perpendicular momentum is relativistic. The resonant absorption condition is 


B' = 7^0 (1 - /3/io) 


(16) 


where 7 is the particle Lorentz factor, Eq is the lab frame energy (in units of mc^) of the radio 
photon, /3 = (1 — 1 / 7 ^)^/^, B' = B/Bcr is the local magnetic held in units of the critical held 
strength B^ = 4.4 x 10^^ G, /xq = cos 6 ^ 0 , and 6q is the angle between the photon direction 
and the particle momentum in the lab frame . In pulsar magnetospheres, particles having 
large Lorentz factors will see radio photons at the cyclotron resonance at high altitude above 
the neutron star surface, when the local magnetic held has dropped low enough to satisfy 
the resonant condition, (from equation [l 6 ]). 


iR 


= 2.8 X 10® 


Bs 

^0,GHz(l — /3ho) 


( 17 ) 


where /xq is the incident absorption angle. 


The particle initially undergoes absorption in low Landau states, where the cyclotron 
emission rate is well below the absorption rate. Therefore, the particle Landau state (and 
pitch angle) will increase stochastically but continuously until the increase in pitch angle 
through resonant absorption and the decrease in pitch angle by synchrotron emission reaches 
an equilibrium. This balance is not reached until the particle occupies a high Landau state, 
where it will radiate synchrotron emission. 


We use the method of HSDF08, based on the work of LP98, to simulate the resonant 
cyclotron absorption and subsequent synchrotron emission from both primary electrons and 
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pairs. LP98 identified two regimes of particle pitch angle increase in the resonant absorption. 
In the hrst regime, when ip 9 q (the particle pitch angle, ip, is less than the radio photon 
incident angle of the, the pitch angle is increasing while the momentnm is nearly constant. 
In the second regime, when {Oq — 'iP) ^ Oq, the pitch angle is constant as the total momentnm 
increases. In the regime where ip <^9 q (eqnation (2.17) of Petrova (2002)), the mean square 
of the pitch angle increases as. 


/ ao{r]')dr]P 

driR 


(18) 


where rj = r/R, rj^ is the radio emission altitude and 

27i^e'^{l -/3fio)Io f eoj{l -/3fio)Y 

-(- R - ) ’ 


(19) 


Here Jq is the observed intensity of the radio emission, in erg ■ cm~‘^ ■ ■ Hz~^ and v is 

its spectral index. Thus, the change in perpendicular momentum due to cyclotron resonant 
absorption is 


dp± 

dt 


abs 


o 1 , (dp 

= 2 ao[g) c— + —\ — 
P± p \dt 


abs 


( 20 ) 


where we used the assumption that p± is proportional to the root mean-square of the pitch 
angle, p± = pippY^"^■ We also compute the evolution of pYY^‘^ rather than computing the 
evolution of the full particle distribution function. 


Combining Eqn (19) and (20), the resonant absorption rate is 


dpi. 

dt 


abs 


T 

= /i—+ 




Pi. 7 


dp 

dt 


abs 


1 <lR 


where 


D = 5.7 X lO^s-^R^ 


^kpc 


<I)o[niJy] (1 - Pdpo), 


( 21 ) 


( 22 ) 


and we have neglected the {dp/dt^^ term in Eqn (21) since the (dp/dt) from acceleration 
and from curvature and synchrotron losses are much larger. In Eqn ( [2^ , we take Jq = 
^o^radd"^/A, where Qrad ~ A/r"^ is the solid angle of radio emission, with A the cross- 
sectional area and r the radius at absorption, <I)o is the measured radio flux (in mJy), and 
d is the source distance (in kpc). If and when the particles satisfy the resonance condition. 


p < Pr (where pr was dehned in Eqn (17)), as they advance along their trajectories, the 


resonant absorption term Eqn (21) will turn on 
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In the regime where Oq — ifj 6 q, we have assumed that the pitch angle maintains a 
constant value of ip = 6 q/2 and that the mean of the particle total momentum is 


P = 



(3 — u)ai 




( 23 ) 


(Petrova 2003), where 


ai (h) 


47r2eV'i(l)Jo 


B'mc J ’ 


V > Vr, 


( 24 ) 


bj = 2e'‘B'ysK‘c and J{ is the derivative of the Bessel function. 

The radio emission is modeled as an empirical cone and core beam centered on the 
magnetic pole to determine the flux distribution of radio photons S{9,eR), given by Eqn (9) 
of HSDF08, where 6 is the magnetic polar angle and Er is the radio photon energy in units 
of mc^. We assume this core component of the flux is emitted at a radius of 1.8R and the 
conal component at radius (in units of the stellar radius) 

0.07 

(25) 


tkg ~ 40 


P 


10 


-15c 


(Kijak & Gil 2003), where Eghz is the radio frequency in GHz. In order to evaluate the radio 
photon intensity <I)o and the incident absorption angle /xo at the position of the particle, 
needed in Eqns (22) and (24), we use this radio core/cone beam model. To evaluate <l)o, we 
divide the radio emission beam in ovc coordinates, dehned in into “beamlets” centered on 
tangents to the held lines and having opening angle ^bm- One set of beamlets is modulated 
by the cone beam at radius tkg (see Eqn (25)) and another set of beamlets is modulated 
by the core beam at radius 1.8R. The absorption angle /xq between each beamlet and a 
particle with which it will interact is determined by computing the travel time of the radio 
photons from their emission points {xbm,ybm,Zbm) to the particle position {xp,yp,Zp), taking 
into account the rotation of the held line during the photon transit time. A particle at a 
given location in the outer magnetosphere can absorb only a fraction of beamlet photons. 
Gontributions to the cyclotron absorption from all the separate beamlets is summed at each 
particle position. 


We have used this empirical model of radio core and cone emission to model the cyclotron 
resonant absorption, although some of the pulsars we study show phase alignment of their 
main radio and high-energy pulses, indicating that the radio emission comes from high 
altitude near the y-ray emission. However, the radio cone emission altitude from the Kijak 
& Gil (2003) model for the Grab at 400 MHz is about 0.27 Plc and for B1937-I-21 is about 
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0.55 Rhc^ which is at fairly high altitude and would likely form a caustic emission pattern to 
give the phase alignment with the gamma-ray pulses. Harding et ah (2008) modeled both a 
cone beam at this single altitude and a cone beam with some extension along the field lines 
to higher altitude for the Crab in and the resulting SR was similar for both cases. In the case 
of the Crab and some MSPs like B1821-I-24 and B1957-I-20 some radio peaks are in phase 
with the gamma-ray peaks and some are near the phase of one of the magnetic poles (as 
for the Crab precursor), indicating that there is both low altitude cone or core emission and 
high altitude radio emission simultaneously. The cones may also be elliptical, not circular, 
from studies of processing pulsars (Weisberg & Taylor 2002). We believe our radio model 
is sufficient for the present calculations since it does place the cone beam emission at high 
altitude for the Crab and the MSPs, in the same region with the gamma-ray emission. 

Although the electron-positron pairs from the PC cascades are thought to produce the 
observed radio emission, we have neglected any energy loss of the pairs due to emission of 
the radio beam. This is reasonable since the radio emission power is a small fraction of 
the spin-down power (about < 10 “®) compared to the fraction of spin-down power used to 
produce the pairs (about 10“^, HMll). 

At each step along its trajectory, the particle radiates an instantaneous synchrotron 
spectrum (Tademaru 1973), 

2^/3 

NsR{e) = ^aR'sinV'£“2/^£;y3exp(-£/£^^), (26) 

where sin-^ = p±/p, p^ = 7 ^ — 1 and £ 3 ^ = ( 3 / 2 ) 7 ^ R'sinis the synchrotron critical 
frequency. 


4.3. Synchrotron Self-Compton Radiation 


Calculation of the SSC emission, the inverse-Compton radiation of both pairs and pri¬ 
mary particles scattering the SR from both pairs and primaries, is carried out in two stages. 
Since the SR photon density is needed to compute the SSC emission, the first stage cal¬ 
culates the SR from both pairs and primary particles at each step along all trajectories in 


the open field volume according to 14.2, and also accumulates the SR emissivity at each 
location, rem = {xem,yem, Zem), In Cartesian coordinates, and photon emission direction, 
? 7 em = {,Vem,x,Vem,y,Vem,z) In the lOF. At each position along a particle trajectory, depend¬ 
ing on whether the particle is a primary electron [p) or a pair (-I-), the SR emissivity is 
incremented by 

Nsnie) Ai 


1 ^em: ^em) 


c Ax ^YtiAy ^ytiAz 

em 


(27) 





15 


where A£ is the spatial step along the particle trajectory in the lOF, rip is the flux of primary 
electrons, given by Eqn Q, is the flux of pairs from Eqn (|^ summed over pair energies, 
and Axem, Ayem, Azem are the spatial divisions in the rem array. At the end of this first 
stage, the SR emission and the total SR emissivity, egji{e,rem,'nem), throughout the open 
magnetosphere of one hemisphere (northern) has been computed. We need the SR emissivity 
from both hemispheres, since particles along trajectories in one hemisphere can scatter SR 
photons produced by particles on trajectories in the other hemisphere, the emissivity from the 
southern hemisphere is computed from the one in the northern hemisphere using reflection 
symmetry: 

I'eirn ^em) I'emj ^em)- (2^) 

The total SR emissivity is then, esnie, rj^rn) = rem, riem) + rem, riem)- 


In the second stage of the simulation, the particles follow their trajectories for a second 
time, radiating SSC using the SR emissivity computed in the first stage. The scattered 
photon distribution from a single particle is. 


dNsscj^s) 

dSsdt 


J dl} J denphie.l^) (1 - P/i) 


(29) 


where nph{e,fl) is the SR photon density, cT(e:,n) is the scattering cross section, and £ and 
Eg are the incident and scattered photon energies. At each step along the particle trajectory, 
the SR photon density in all directions is computed at the current position, fiof, of the 
particle. 


rioF, r/em) ^ / dXg 


rymax 

' ymin 


dye 


dZe 


^em; 'He 


where 

= {XlOF — XemY + {VlOF — UemY + {^lOF — ^erriT ■ 


(30) 


(31) 


In principle, the Klein-Nishina (KN) scattering cross section should be used in Eqn 
(29). In practice, a full integration over all angles using the full KN cross section at each 


particle position is not computationally feasible. We have therefore adopted an approximate 
approach, using the formula of Jones (1968) for the KN photon production rate for a sin¬ 
gle ultra-relativistic particle with Lorentz factor 7 scattering an isotropic distribution, and 
modulated by the anisotropic photon density and the relative velocity factor: 


d£ssci^ ^ 1 f 

de.dt 4 f 


de nph{e,rioY,r] em ) (1 - 


dnKN{£,es) 

de.de 


(32) 


where 
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dnKN{e,es) 27rrgc , n i o ^ i ^ M 

-= [2glng+(l + 2g)(l-g) +- ^ ^ ^ J l - g)], 


de.de 




21 + rg 


(33) 


r = 4 £ 7 , q = 


El eg 

- Et = — 

m-EiY ' 7 ’ 


(34) 


and ro is the classical electron radins, fi = cos^ is the incident photon angle and (3 = 


(1 — 1 / 7 ^)^/^. In the limit F << 1 and El « 1, Eqn ( |32[ ) rednces to the Thompson limit 
resnlt (Blnmenthal & Gould 1970). 

To compute the SSC spectral contribution from a particle at position fiof with velocity 


V in the lOF, the integration in Eqn (32) is performed by looping over the whole range 


of incident photon energy e and directions p and 0, relative to the particle direction. We 
assume that the scattered photon direction, fig and Ys, is the same as the particle direction, 
r/e, so 


, , I Ve,y 

l^s = Ve^z, tan 0 s = - 


(35) 


To evaluate the photon density for each incident photon direction we need to know the 
incident photon direction, /ij and 0*, in the lOF, 


/ij = psh + sin 9g sin 9 cos 0 


and 


with 


cosA 0 j = 


04 0 s “t“ A 02 

COS 6* — cos 6*0 cos 6*,- 


(36) 

(37) 

(38) 


sin 9g sin 9i 

The quantities p, and 0j dehne the direction 77 * in the lOF that is then used to hnd the SR 
photon density in that direction. 


4.4. Particle Dynamics 


The equations of motion for the Lorentz factor, 7 , and perpendicular momentum, p± 
(in units of me) of a particle as it moves along a held line can be written (Harding, Usov & 
Muslimov 2005),f 

dt me 3m^e^ 3pl \dt / \dt / 
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^^_3c 2 e^ q 2P± f dp±{l) \ 

dt 2 3m^c® 7 \ dt J 


(40) 


In equation (39), the different terms of the right hand side are acceleration, synchrotron 
losses, curvature radiation losses, cyclotron/synchrotron absorption and inverse Compton 
losses. In equation (40), the various terms of the right hand side are adiabatic momentum 
change along the field line, synchrotron losses and cyclotron/synchrotron resonant absorp¬ 
tion. The SSC losses are negligible for p±. 


In the first stage, integrating the particles trajectories to compute the SR and the SR 
emissivity, only SR losses are included since the SSC losses are computed in the second 
stage. However, we find that the SSC losses are much smaller than the SR losses in all cases 
except for the very high altitude part of the pair trajectories for the Vela pulsar (see Figure 
[^discussed below). For this case, discussed in the SSC losses are negligible, and the SR 
losses do not change the particle energy significantly since the cyclotron absorption balances 
the SR losses. The spatial step size for the particle trajectories is dynamically adjusted at 
each step to not exceed limits set by radius of curvature of the trajectory (to limit fractional 
change in emission direction to 10 %), acceleration rate (to limit fractional gain in primary 
Lorentz factor to 10%), CR and SR loss rates (to limit fractional energy loss to 10%), and 
absorption rate (to limit fractional increase in p± to 50%). 


5. Results 

The radiated spectra of primary particles and pairs have been simulated for the Crab 
pulsar, the Vela pulsar, and two of the most energetic MSPs in the Northern and Southern 
hemispheres, PSR B1937-I-21 and B1821-24. The pair spectra computed for each of the 
pulsars using the cascade simulation described in are used to compute the spectra of SR, 
CR and SSC radiation. However, we tried several different pair multiplicities as well as the 
multiplicity from the original cascade simulation. 

The particle trajectories start at the neutron star surface and extend to rmax = 2.0i?LC 
in radius or = 1.9i?LC in cylindrical radius, whichever is reached first. The limit 

is set so that the calculation stays inside the grid of the numerical field lines. The region of 
footpoints of injected particles on the PCs is divided into 5 rings between and each 
of which is divided into 15 azimuthal divisions (Aazim = 0.04 divisions per degree), totaling 
75 each of injected primary particles and pairs at each of 34 energies. The code was run on 
75 parallel processors on the Discover cluster at Goddard. For the calculation of the SSC 
spectrum, we used 20 divisions in cos d and 12 divisions in cj) in the integration over incident 
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photon directions (Eqn [^). The SR photon emissivity array has 50 logrithmic energy 
bins, from 2 x 10“® — 10^ MeV, 18 angular divisions in each of three independent Cartesian 
direction cosines, r/em = {Vem,x,Vem,y,Vem,z) from —1 to 1, and 9 divisions in each of three 
spatial coordinates rem = {xem,yem, Zem) from — 2.0i?LC to 2.0-Rlc- To assure that there 
was enough angular accuracy in the SR photon emissivity array of the SSC computation, we 
tested the computed SSC flux level with increasing number of angular divisions. It was found 
that the SSC flux decreases with increasing number of divisions but stabilizes at 18 divisions 
in each component of r]em, for a combination of 18^ = 5832 elements over dvr steradians, so 
this number was adopted. 


Figure shows examples of the evolution of the pair perpendicular momentum, p±, 
cyclotron absorption rate, {dp±/diY^^, synchrotron loss rate, {d'y/d£)^^, SSC loss rate, 
(dj/dl)^^^, and synchrotron critical energy, esr, as a function of radial distance along a 
pair trajectory for the Crab, Vela and B1937+21, from integration of Eqns (39) and (40), all 
for pair energies around 7 + ~ 10^. The cyclotron absorption starts at the radius of the radio 
cone emission, which is about 0.2i?LC for fhe Crab, about O.IRlc for Vela, and 0.56 Rlc for 
B1937+21, according to Eqn (25). Before the radio emission altitude is reached, there is 
some transient SR from having to set a small non-zero initial p± = 10“®, to avoid overflow 
in evaluation of {dp±/d£)°'’^^ in Eqn (21). The fluctuations in the absorption rate as a func¬ 
tion of radius are due to emission from different beamlets the electron encounters (see ^4.2), 
reflecting the numerical resolution of these sub-elements of the radio cone beam. For this 
pair energy, the p± becomes relativistic. In the case of the Crab and B1937-I-21, the radio 
photons are emitted at high enough altitude and the magnetic field strength near the light 
cylinder, Rlc, is high enough that the radio photons stay in resonance in the particle rest 
frame through the entire particle trajectory. Their Lorentz factors (not plotted) are nearly 
constant since the cyclotron absorption rate comes into equilibrium with the SR loss rate. 
In the case of Vela, the absorption occurs over only a small part of the trajectory above 
the radio emission altitude, above which the particle is no longer “in sight” of the relatively 
more narrow radio beam. The p±, {d'j/d£)^^ and Esr thus drop with increasing radius. The 
SSC loss rate is always less than the SR loss rate, except in the case of Vela where the SR 
loss rate drops below the SSC loss rate at higher altitude. The SSC loss rate for Vela does 
not decrease much at high altitude. Because the Rlc for Vela is much lower compared to 
the other pulsars, the SR critical frequency is about four order of magnitude lower, so that 
the radiated photon number density increases with altitude. Although the SSC loss rates 
are not included in the trajectory integration, they are never large enough to change the 
particle energy significantly. 


Figure shows examples of the evolution of the primary electron Lorentz factor 7 , 
curvature radiation loss rate, {d'y/d£)^^, synchrotron loss rate, {d'y/d£)^^ and SSC loss rate. 
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, as a function of radial distance along a pair trajectory for the Crab, Vela and 
B1937+21. For the primaries, CR losses are always much larger than SR and SSC losses and 
the electrons reach CR reaction limit in all cases within a few tenths of Rhc- The maximnm 
Lorentz factors for the chosen valnes of acceleration rate are slightly over 10^, giving critical 
CR energies of a few GeV. The SR loss rate peaks soon after the start of the cyclotron 
resonant absorption, at which point the particle Lorentz factor drops as the SR loss and 
acceleration gain rates are briefly in eqnilibrinm (Harding et al. 2005). As the SR loss rate 
drops, the particle then retnrns to eqnilibrinm between acceleration gain and CR loss rates. 


5.1. Crab Pulsar 

For the Crab pnlsar, we use the following measnred parameters: P = 0.033 s, P = 
4.22 X 10“^^ss“^, d = 2 kpc and S'400 = 700 mJy, where d is the distance and S'400 is the 
radio flnx at 400 MHz. In addition, we assnmed a magnetic inclination angle of a = 45°, a 
constant acceleration rate of Race = eE\\/mc^ = 0.2 cm“^ = 6x 10“^ Blc, where Rlc = 6x 10® 
G, for the primary particles, and two pair mnltiplicities, M+ = 2 x 10^, from a steady-state 
pair cascade, and M_|_ = 3 x 10®, from a time-dependent pair cascade. 

Fignre shows the observed and modeled spectral energy distribntion (SED) of the 
phase-averaged flux of the Grab pulsar from optical to VHE y-ray energies at viewing angle 
C = 60° and for M+ = 2 x 10^. The different model radiation components, primary GR, 
SR, SSG and pair SR and SSG, are plotted separately. We mnltiplied the phase-averaged 
pair SR flnx by a factor of 11 to match the observed spectral points and mnltiplied all the 
other components by the same factor. The shape of the pair SR component, dependent 
on the shape of the pair spectrnm, and its peak energy, dependent on the magnetic held 
strength near the light cylinder, nicely matches the observed spectrnm from optical throngh 
the BeppoSax hard X-ray measnrements. However, it falls short of the GOMPTEL low- 
energy y-ray points. The pair SSG component peaks at several GeV bnt falls well below the 
observed Fermi and VHE spectral points. The primary GR component peaks near a GeV 
bnt way above the observed spectrnm, while the lower primary SR component peaks aronnd 
1 MeV at the level of the observed spectrnm. The primary SSG component, peaking sharply 
at a few times 10® MeV, has a hnx below the scale of the plot. 

Fignre shows the observed and modeled SED for the same parameters as the model 
in Fignre |^bnt for a higher pair mnltiplicity of M+ = 3 x 10®. For this case, the pair SR 
spectrnm matches the observed spectrnm withont any renormalization factor. The associated 
pair SSG component now ronghly matches the level of the higher-energy Fermi points and 
the MAGIG and VERITAS measnrements. For the adopted valne of the acceleration rate 
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Ra.cc -1 the primary CR component fills in at a few hundred MeV. However, the primary SR 
component is too low to account for the observed soft 7 -ray emission. The tip of the steeply 
rising, nearly mono-energetic primary SSC spectrum is just visible at a few TeV at the 
bottom right edge of the plot, and the pair CR lies several orders of magnitude below the 
plot lower boundary. 

We tested the addition of a high-energy power law extension to the pair spectrum (shown 
in Figure]^ whose SR spectrum would account for the observed emission in the 1 - 100 MeV 
range. The resulting pair SR and SSC opponents are shown in Figure as dashed lines. 
The SSC spectrum now extends to higher energy with a harder spectrum that exceeds the 
observed MAGIC and VERITAS points. This extension of the pair spectrum, which has no 
physical basis, would thus produce an SSC spectrum that seems to be ruled out by the data. 
This implies that the observed 1-10 MeV emission is not produced by the same particles 
that produce the SSC emission. 


5.2. Vela Pulsar 

For the Vela pulsar, we use the following measured parameters: P = 0.089 s, P = 
1.25 X 10“^^ss“^, d = 0.29 kpc and 5400 = 5000 mJy. We assumed a magnetic inclination 
angle of a = 75°, a constant acceleration rate of Race = eE\\/mc^ = 2.0 = 0.08Rlc for 
the primary particles, and two pair multiplicities, M+ = 6 x 10 ^, from the steady-state pair 
cascade, and M+ = 1 x 10 ^, for a time-dependent pair cascade. 

Figure shows the observed and models SED of the phase-averaged flux of Vela at 
viewing angle C, = 60° for the two different values of M+. We multiplied the phase-averaged 
primary CR flux by a factor of 0.14 to match the level of the observed Fermi spectral points 
and multiplied all the other components by the same factor. The value of Race had been 
chosen so that the SED peak of the primary CR spectrum matches the peak of the Fermi 
spectrum, but the flux level being too high implies that either the flux of primary particles 
or the distance over which they are radiating is smaller than what we assume. The pair SR 
and SSC components peak at energies of around 1 keV and 1 GeV respectively, somewhat 
lower than the peaks of those components for the Grab. For both values of M_|_, the SSG 
flux is orders of magnitude lower than the primary GR flux. The primary SSG component 
is well below observable levels. The striking differences between the Vela and Grab SEDs 
come from a number of differences in their intrinsic properties. First, although they have 
similar surface magnetic field strengths, the larger magnetosphere of Vela gives it a much 
lower Rlc = 4 x 10^ G compared to that of the Grab. This produces the lower peak energies 
of the SR and SSG components and also the lower flux levels because the magnetic field 
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drops below the value where the radio frequencies are in the cyclotron resonance in the 
pairs’ rest frames well before the light cylinder, whereas for the Crab, the radio photons are 
in resonance out to and beyond the LC. Second, the radio luminosity of Vela is lower by 
about a factor of 5 compared to the Crab, and is produced at lower altitude relative to the 
LC. The particles therefore see a lower density of radio photons over their trajectories and 
undergo less absorption. Third, the pair spectrum of Vela has a high-energy turnover at a 
lower energy, producing the lower peak energies of SR and SSC spectra. 


5.3. Millisecond Pulsars 

MSPs are promising sources of high non-thermal SR and SSC emission, since they have 
small magnetospheres and consequently a number have large Rlc- In fact, a larger fraction 
of MSPs have observed hard non-thermal X-ray spectra that do non-recycled pulsars. We 
model two of the most energetic MSPs, B1821-24 and B1937-I-21, that have the highest Rlc 
and high levels of non-thermal X-ray emission. B1821-24 (J1824-2452), with P = 3.05 ms, 
P = 1.6 X 10“^®ss“^ and 5*400 = 40 mJy, lies in the globular cluster M28, at a distance of 
d = 5.5 kpc. B1937+21 (J1939+2134), with P = 1.6 ms, P = 1.0 x 10-l9ss-^ 5400 = 240 
mJy and d = 3.6 kpc, was the hrst discovered MSP (Backer et ah 1982). Both are isolated 
MSPs, have Plc = 7-3 x 10® G and 10® G respectively, exhibit giant radio pulses (Bilous et 
al. 2015, Popov & Stappers 2003) and have aligned high-energy and radio prohle components 
similar to the Crab. 

Figure [^shows the observed and modeled spectra of B1821-24 from soft X-rays through 
VHE 7 -ray energies for a = 45° and viewing angle ( = 80°. Although there are no strong 
constraints on either a or (, modeling the y-ray light curve only favors a ~ 40° and ( = 85° 
(Johnson et al. 2013). We compute SR and SSC spectra for M+ = 10^, a value for a pair 
cascade with offset PC of e = 0.6, and also an extreme value of M+ = 10®, both assuming 
Race = eE\\/mc^ = 2.0 = 5 x 10 “^Plc for the primary particles. The pair SR component 
peaks between 1 and 10 MeV and the pair SSC component peaks around 50 GeV, almost 
two decades higher than the SR and SSC components in the Crab SED. This is expected, 
given than the pair spectrum of B1821-24 is higher in energy by a similar factor relative 
to the Crab pair spectrum, and the Plc values are similar. The pair SR spectrum for 
the M_|_ = 10® case matches the slope and level of the X-ray data with a minor shift by a 
factor of 0.8, which has been applied to the other components for this model. The primary 
CR component roughly coincides with the Fermi points, with the primary SR being much 
lower. The pair SSC spectrum lies about four decades below the pair SR and primary CR 
components and is not detectable. The pair SR flux for M+ = 10® is almost two decades 
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above the observed X-ray flux, although the pair SSC flux might be detectable with Fermi 
and even HESS. Since time-dependent cascade simulations for MSPs are not available at 
present, we cannot say whether such high multiplies are in fact achievable, although a pair 
multiplicity much higher than M_|_ = 10^, which already accounts for the observed X-ray 
spectrum, seems to be ruled out. Therefore, for our chosen values of inclination and viewing 
angle, we do not predict a detectable VHE component from B1821-24. 

Figure]^ shows the observed and modeled spectra of B1937-1-21 from soft X-rays through 
VHE 7 -ray energies for a = 75° and viewing angle ( = 70°. Johnson et ah (2014) fit the 
combined y-ray and radio light curves of this MSP using several geometrical emission models. 
For the two-pole caustic geometry, they And a = 88°(-l-2 — 1)°,C = 88°(-|-2 — 1)°, and for 
an OG geometry, a = 72° ± 1°,(C = 85° ± 1°. We show the spectra for pair multiplicities 
of M+ = 10^ of a pair cascade with e = 0.6, and also for M_|_ = 10^, both assuming 
7?acc = 2.0 = 3 X 10“^ Hlc for the primary particles. Applying only the l/(27rd^) factor to 
obtain the phase-averaged flux from the sky map, the pair SR flux for M_|_ = 10® is consistent 
with the Chandra X-ray flux. The peak of the pair SR SED is around 5 MeV, very similar 
to the SR SED of B1821-24. The primary CR component is about a factor of 3 below the 
Fermi flux points, and the pair SSC component, peaking at around 100 GeV, is about a 
factor of 2 below CTA sensitivity. The pair SR component for the case of M+ = 10® is about 
two decades below the observed X-ray flux and the pair SSC emission lies well below the 
plot boundary. The SR fluxes of B1821-24 and B1937-I-21 for the same pair multiplicity and 
similar are very different since the model inclination angles are quite different. The value of 
a = 75° was chosen to match the best model fits of the y-ray light curve, but a simulation 
for B1937-I-21 with a = 45° would match the observed X-ray spectrum with a lower pair 
multiplicity. In this case, the SSC component would be much lower. 


6. Discussion 

We have simulated the synchrotron and synchrotron-self Compton emission from a broad 
spectrum of electron-positron pairs in a global force-free magnetosphere. The scattering takes 
place in the extreme Klein-Nishina limit and although the angular dependence of the cross 
section is neglected, the full 3D angular dependence of the SR photon density is treated by 
storing the SR emissivity over the entire magnetosphere from both rotational hemispheres. 

Our simulation of the Crab pulsar radiation reproduces both the flux level and the 
shape of the observed optical to hard X-ray emission assuming a pair multiplicity of M+ = 
3 X 10®. Such a high multiplicity is about an order of magnitude larger than that produced in 
steady-state PC cascades, but is achievable with the more realistic time-dependent cascades 


(Timokhin & Harding 2015). The predicted SSC flux in this case roughly matches the tail 
of the observed Fermi spectrum as well as the MAGIC and VERITAS detected points. A 
CR component from primary particles is necessary to explain the Fermi spectrum below 
a few GeV. The model SSC spectrum is not a power law, but reflects the shape of the 
pair spectrum, which has a smooth and continuous curvature extending to ~ 1 TeV. An 
SSC component from primary particles appears above a few TeV and may ultimately be 
detectable. 

The Vela pulsar does not produce detectable SSC emission in our simulations, even for 
a pair multiplicity of 10®. The predicted SSC emission has both a lower flux and a lower 
SED peak energy than the Crab SSC emission, because the SED peak energy of the Vela 
pair SR is several decades lower. This difference is due to a combination of lower pair energy 
and lower Rlc- However, the Fermi spectrum can be produced by primary CR and SR. 

Simulating emission from two of the most energetic MSPs, B1821-24 and B1937+21, we 
And that their observed X-ray emission is well reproduced by SR from pairs. The SEDs of 
MSP SR components are predicted to peak at 1 - 10 MeV and could be detected by proposed 
Compton telescopes that would test this model. Since the soft X-ray emission is produced 
by the pairs at the low-energy turnover in the pair spectrum, well below the SED peak, the 
spectral index is that expected for SR from a single particle, —2/3 (4/3 for the SED). This 
index is much higher than the soft X-ray index of the pair SR spectrum of the Crab, where 
the soft X-ray emission is near the SED peak and produced by a range of pair energies. 
However, the predicted SSC fluxes of the MSPs are too low to be detectable with current 
instruments. 

Thus from our simulations, only Crab and Crab-like pulsars, such as B1509-58 and 
B0540-69 in the LMC, are expected to have high enough levels of SSC emission at VHE 
energies to possibly be detectable by ground-based Air-Cherenkov telescopes. They share the 
characteristics that are important for production of strong SSC emission: high levels of non- 
thermal X rays, high pair multiplicity and high magnetic fields near the light cylinder. Non- 
Crab-like, middle-aged pulsars like Vela have much lower levels of non-thermal X rays relative 
to their GeV emission. So even if they are able to generate very high pair multiplicities ~ 10®, 
their pair energies will be too low and their SSC emission will not reach detectable levels 
at VHE energy. MSPs have higher levels of SSC than middle-aged pulsars, and their SSC 
spectra extend to higher energies than Crab-like pulsars, since their pair spectra reach TeV 
energies. However, their pair multiplicities may not be high enough to generate fluxes of 
SSC emission that are detectable by current telescopes. These findings are consistent with 
the recent stacking analysis of McCann (2015). 

Lyutikov (2013) presented a model for the Crab pulsar emission, assuming cyclotron 



and cyclotron-self Compton emission by pairs prodnced in the OG and gyrating at the 
cyclotron resonance with small pitch angles. It was snggested that the pairs acqnire their 
pitch angles throngh coherent emission of radio waves, bnt their perpendicnlar momenta 
are assnmed to remain non-relativistic so that the radiation occnrs in the cyclotron rather 
than the synchrotron regime. The reqnired pair mnltiplicity is M+ = 10® — 10^, well above 
what can be prodnced in either PC or OG pair cascades. In onr model, the pairs gain pitch 
angles throngh resonant absorption and acqnire relativistic perpendicnlar momenta, so that 
the radiation is in the synchrotron regime. Since the synchrotron radiation power is mnch 
higher than cyclotron power, we can prodnce the Crab optical to X-ray emission with a mnch 
lower pair mnltiplicity of 3 x 10®. 
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Fig. 1.— Pair spectra (in units of pairs/s/mc^) from polar cap pair cascade simulations for 
the Crab, Vela, B1821-24 and B1937+21 pulsars, as labeled. The dashed line is a power law 
extension to the Crab pair spectrum (see text). 
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Fig. 2.— Evolution of the dynamics and radiation of pairs at initial energy 7 + as a func¬ 
tion of radius (in units of light cylinder radius) along their trajectory, for the Crab, Vela 
and B1937-I-21. Quantities plotted are perpendicular momentum, p± (in units of me), syn¬ 
chrotron loss rate, (cm“^), SSC loss rate, (dy/df')®®'" (cm“^), cyclotron absorption 

rate, {dp±/d£)’^^^ (cm“^), and critical synchrotron energy £ 3 ^,^ (in units of me?). 
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Fig. 3.— Evolution of the dynamics and radiation of primary electrons as a function of radius 
(in units of light cylinder radius) along their trajectory, for the Crab, Vela and B1937+21. 
Quantities plotted are Lorentz factor, 7 , synchrotron loss rate, (cm“^), CR loss 

rate, {dj/di)^^ (cm“Q, SSC loss rate, (dj/di)^^^ (cm“^), and acceleration gain rate, Race 
(cm“^). 
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Fig. 4.— Model spectra of phase-averaged pulsed emission components from primary elec¬ 
trons and pairs (as labeled) from the Crab pulsar, for magnetic inclination angle a = 45° 
and observer angle C, = 60° and pair multiplicity = 2 x 10"^. Data points are from Kuiper 
et ah (2001) |http://www.sron.nl/divisions/hea/kuiper/data.html|, Abdo et al. (2013) 
[http://fermi.gsfc.nasa.gov/ssc/data/access/lat/2nd_PSR_catalog/], Aleksic et al. (2012) 
and Aliu et al. (2011). 
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Fig. 5.— Model spectra of phase-averaged pulsed emission components from primary elec¬ 
trons and pairs (as labeled) from the Crab pulsar, for magnetic inclination angle a = 45° 
and observer angle C, = 60° and pair multiplicity M+ = 3 x 10^. The dashed lines are the 
SR and SSC spectra resulting from a power law extension to the cascade pair spectrum. 
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Energy (GeV) 

Fig. 6.— Model spectra of phase-averaged pulsed emission components from primary 
electrons and pairs (as labeled) from the Vela pulsar, for magnetic inclination angle 
a = 75° and observer angle C, = 60°. Solid lines are for pair multiplicity M_|_ = 
6 X 10^ and dashed lines are for M+ = 10^. Data points are from Abdo et ah 
(2013) |http://fermi.gsfc.nasa.gov/ssc/data/access/lat/2nd_PSR_catalog/| and Harding et 
al. (2002). 
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Fig. 7.— Model spectra of phase-averaged pulsed emission components from primary elec¬ 
trons and pairs (as labeled) from PSR B1821-24, for magnetic inclination angle a = 45° and 
observer angle ( = 80°. Solid lines are for pair multiplicity M_|_ = 10^ and dashed lines are 
for M_|_ = 10^. Data points are from Kuiper et al. (2005) and Johnson et al. (2014). The 
thick pink and yellow lines are the HESS and CTA sensitivity limits, respectively. 
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PSR B1937+21 



Energy (MeV) 


Fig. 8.— Model spectra of phase-averaged pulsed emission components from primary elec¬ 
trons and pairs (as labeled) from PSR B1937+21, for magnetic inclination angle a = 75° 
and observer angle ( = 70°. Solid lines are for pair multiplicity M+ = 10^ and dashed lines 
are for M+ = 10^. Data points are from Ng et ah (2014) and Guillemot et ah (2012). The 
thick yellow line is the CTA sensitivity limit. 













